#---- Supplementary replication script for  ----

#Turnbull-Dugarte, López-Ortega & Hunklinger (2024)
#"Do citizens stereotype Muslims as the Illiberal 'Bogeyman'? Evidence from a Double-List Experiment"

##This file produces the summary statistics in the appendix: Tables A1-A4##

remove(list=ls())

#install.packages("pacman")
library(pacman)
p_load(here, tidyverse, readr, modelsummary, estimatr, haven)

set.seed(123)

list_dataUK <- read_dta("data/listUK.dta")
list_dataUK<- list_dataUK  %>% 
  mutate(treated=as.numeric(treated),
         value=as.numeric(value),
         respondent_gender=as.numeric(respondent_gender))
list_dataUK_singleround <- subset(list_dataUK, round == 1)

list_dataUS <- read_dta("data/listUS.dta")
list_dataUS<- list_dataUS  %>% 
  mutate(treated=as.numeric(treated),
         value=as.numeric(value),
         respondent_gender=as.numeric(respondent_gender))
list_dataUS_singleround <- subset(list_dataUS, round == 1)


list_dataDE <- read_dta("data/listDE.dta")
list_dataDE <- subset(list_dataDE, !(is.na(weight)))
list_dataDE<- list_dataDE  %>% 
  mutate(treated=as.numeric(treated),
         value=as.numeric(value),
         respondent_gender=as.numeric(respondent_gender))
list_dataDE<- as.data.frame(list_dataDE)
list_dataDE_singleround <- subset(list_dataDE, round == 1)


list_dataNL <- read_dta("data/listNL.dta")
list_dataNL <- subset(list_dataNL, !(is.na(weight)))
list_dataNL<- list_dataNL  %>% 
  mutate(treated=as.numeric(treated),
         value=as.numeric(value),
         respondent_gender=as.numeric(respondent_gender))
list_dataNL<- as.data.frame(list_dataNL)|> drop_na(weight)
list_dataNL_singleround <- subset(list_dataNL, round == 2)



###SUMMARY STATS###

datasummary(treated + respondent_age + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality 
            ~ N + Mean + SD + Min + Max,
            data = list_dataUS_singleround, output="latex")


datasummary(treated + respondent_age + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality 
            ~ N + Mean + SD + Min + Max,
            data = list_dataUK_singleround, output="latex")

datasummary(value + treated + respondent_age + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality 
            ~ N  + weighted.mean * Arguments(w = weight,na.rm = TRUE)+ SD + Min + Max,
            data = list_dataDE_singleround, output="latex")

datasummary(value + treated + respondent_age + respondent_gender + respondent_degree +
              respondent_ideology + affective_muslims + affective_lgb + respondent_sexuality 
            ~ N + weighted.mean * Arguments(w = weight,na.rm = TRUE)+ SD + Min + Max,
            data = list_dataNL_singleround, output="latex")


desired_columns <- c('treated', 'respondent_gender', 'respondent_age', 'respondent_sexuality', 'respondent_degree', 'respondent_ideology', 'affective_lgb', 'affective_muslims')
desired_columnsW <- c('weight', 'treated', 'respondent_gender', 'respondent_age', 'respondent_sexuality', 'respondent_degree', 'respondent_ideology', 'affective_lgb', 'affective_muslims')

summarylistUS_A <- subset(list_dataUS_singleround, select = desired_columns)
summarylistUK_A <- subset(list_dataUK_singleround, select = desired_columns)
summarylistDE_A <- subset(list_dataDE_singleround, select = desired_columnsW)
summarylistNL_A <- subset(list_dataNL_singleround, select = desired_columnsW)

summarylistNL_A <- summarylistNL_A |> 
  mutate(respondent_degree= as.double(respondent_degree))

summarylistDE_A <- summarylistDE_A |> 
  mutate(respondent_degree= as.double(respondent_degree))

datasummary_balance(~treated,
                    data =summarylistUS_A,
                    stars=TRUE,
                    output="latex")

datasummary_balance(~treated,
                    data =summarylistUK_A,
                    stars=TRUE,
                    output="latex")

datasummary_balance(~treated,
                    data =summarylistNL_A,
                    stars=TRUE,
                    output="latex")

datasummary_balance(~treated,
                    data =summarylistDE_A,
                    stars=TRUE,
                    output="latex")

